Crack Path Prediction in Anisotropic Brittle Materials 
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A force balance condition to predict quasistatic crack paths in anisotropic brittle materials is 
derived from an analysis of diffuse interface continuum models that describe both short-scale failure 
inside a microscopic process zone and macroscopic linear elasticity. The derivation exploits the 
gradient dynamics and translation symmetry properties of this class of models to define a generalized 
energy-momentum tensor whose integral around an arbitrary closed path enclosing the crack tip 
yields all forces acting on this tip, including Eshelby's configurational forces, cohesive forces, and 
dissipative forces. This condition is validated quantitatively by numerical simulations. 
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The prediction of the path chosen by a crack as it prop- 
agates into a brittle material has been a long standing 
problem in fracture mechanics. This question has been 
addressed primarily in a theoretical framework where the 
equations of linear elasticity are solved with zero traction 
boundary conditions on crack surfaces that extend to a 
sharp tip The stress distributions near the tip have 
universal divergent forms 
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where Km are the stress intensity factors for the three 
standard modes I, II, or III of fracture (m = 1, 2 or 3) and 
Q is the angle between the radial vector of magnitude r 
with origin at the crack tip and the local crack direction. 
For the crack to propagate, the energy release rate (or 
crack extension force) 



G = a{Kl+Kl)+Kl|{2^x), 
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must exceed some material dependent threshold Gc that 
is theoretically equal to twice the surface energy (Gc = 
27), but often larger in practice. Here, v is Poisson's 
ratio, E is the bulk modulus, ^ is the shear modulus, 
and a = (1 - u'^)/E. 

Like other problems in fracture, the prediction of crack 
paths was first examined 2] for mode III which is sim- 
pler because the antiplane component of the displace- 
ment vector 1*3 is a purely scalar Laplacian field. The 
stress distribution near the tip, can be expanded as 
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and the dominant divergent contribution is always sym- 
metrical about the crack direction. This implies that 
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knowledge of K3 cannot predict any other path than 
a straight one. To avoid this impasse, Barenblatt and 
Cherepanov retained the subdominant term ~ sin0, 
which breaks this symmetry, and hypothesized that a 
curvilinear crack propagates along a direction where 
A2 — 0, and hence when the stress distribution is sym- 
metrical about the crack direction. In subsequent exten- 
sions of this work, several criteria have been proposed for 
plane loading, for which the tensorial nature of the stress 
fields makes it possible to predict non-trivial crack paths 
urely from the knowledge of the stress-intensity factors 
Q. The generally-accepted condition "if2 = 0" as- 
sumes that the crack propagates in a pure opening mode 
with a symmetrical stress distribution about its local axis 
and is the direct analog for plane strain (7/3 = 0) of 
the condition A2 — for mode III. This "principle of lo- 
cal symmetry" has been rationalized using plausible ar- 
guments but cannot be derived without an explicit 
description of the process zone, where elastic strain en- 
ergy is both dissipated and transformed nonlinearly into 
new fracture surfaces. As a result, how to extend this 
principle to anisotropic materials, where symmetry con- 
siderations have no obvious generalization, is not clear 
5]. In addition, path prediction remains largely unex- 
plored for mode III even for isotropic materials. 

In this letter, we address the problem of path predic- 
tion in the context of continuum models of brittle fracture 
that describe both short scale failure and macroscopic 
linear elasticity within a self-consistent set of equations. 
Such models have already proven capable to reproduce 
a wide range of phenomena for both antiplane Q and 
plane loading from the onset of crack propagation at 
the Grifhth threshold to dynamical branching instabili- 
ties and oscillatory instabilities. From an analysis 
of these models, we derive a new condition to predict 
crack paths that is interpreted physically in the context 
of previous results from the fracture community. 

For clarity of exposition, we base our derivation on 
the phase-field approach of Ref. where the displace- 
ment field is coupled to a single scalar order parame- 
ter or "phase field" 0, which describes a smooth tran- 
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sition in space between unbroken = 1) and broken 
states (0 = 0) of the material. Our approach is suffi- 
ciently general, however, to be applicable to a large class 
of diffuse interface descriptions of brittle fracture. We fo- 
cus on quasi-static fracture in a macroscopically isotropic 
elastic medium with negligible inertial effects. Material 
anisotropy is simply included by making the surface en- 
ergy, 7(6*), dependent on the orientation 6 of the crack 
direction with respect to some underlying crystal axis. 

For brevity of notation, we define the four-dimensional 
vector field ijj^ = Uk for 1 < fc < 3 and tjj^ = (j) where Uk 
are the components of the standard displacement field. 
The energy density £ depends on (j) and djip'' = dip'' /dxj, 
where spatial gradients of the displacement contribute to 
the elastic strain energy and gradients of the phase-field 
contribute to the surface energy |^. The equations of 
motion are derived variationally from the spatial integral 
of £, the total energy E of the system, and obey the 
gradient dynamics 
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where 6i^i = 1 and 5ij = ior i j. These Euler- 
Lagrange equations for ip'' = Uk are simply the static 
equilibrium conditions that the sum of all forces on any 
material element vanish. The fourth equation for ip* = cp 
is the standard Ginzburg-Landau form that governs the 
phase-field evolution, where % is a kinetic coefficient that 
controls the rate of energy dissipation in the process zone. 

In the present model that describes both microscopic 
and macroscopic scales, the problem of predicting the 
macroscopic path of a crack can be posed as an "inner- 
outer" matching problem. We seek inner solutions of 
the equation of motion Q on the scale ^ of the process 
zone subject to far field boundary conditions imposed 
by matching these solutions to the standard solutions of 
linear elasticity on the outer scale of the system W ^ 
^. These outer solutions change slowly on a scale where 
the crack advances by a distance ~ ^. We can therefore 
search for inner solutions by rewriting Eq. Q in a frame 
translating uniformly at the instantaneous crack speed V 
parallel to the crack direction 



Sk,4Vx ^di(p = dj 



d£ 



d£ 



ddjip'' dip'' 



(5) 



We then seek solutions of Eq. ^ with far displacement 
fields (r ^ ^) in the unbroken solid that yield the singular 
stress distributions defined by Eqs. ^ and ||2J). 

There are two ways to proceed to solve this problem. 
The first approach, which will be exposed in more de- 
tails elsewhere, exploits the existence of the stationary 
semi- infinite crack solution, ip^, for K2 — A2 — and 
G = Gc which is symmetrical about the crack axis. One 
can then seek solutions of Eq. jS)) linearized around this 
Griffith crack with the driving force for crack advance 



G — Gc and symmetry breaking perturbations (^^2, ^2, 
and the anisotropy of the surface energy dg^) assumed 
to be small. Owing to the variational character of the 
phase-field equations, the linear operator of the resulting 
linearized problem is self-adjoint and hence has two zero 
modes, diip'^, associated with translations of the Griffith 
crack. The standard requirement that non-trivial solu- 
tions to this problem be orthogonal to the null space of 
the adjoint linear operator yields two independent solv- 
ability conditions (for i = 1,2). 

The second approach, which we adopt here, exploits 
the variational character of the equations of motion. It 
yields identical solvability conditions as the first approach 
when G — Gc and symmetry breaking perturbations are 
small, but it is more general since it does not require 
these quantities to be small. We start from the equality 
obtained simply from chain rule differentiation. 
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where dx = dxidx2 and is an arbitrary region of the 
(xi,X2) plane. Using Eq. Q to eliminate d£/dipk from 
the integrand of the right-hand-side (r.h.s.), we obtain 

F,= I dxdjTij-— [ dxdi(pd,(p = for i = l,2. (7) 
Jn X J 



The generalized energy-momentum (GEM) tensor 
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extends the classical energy-momentum tensor of lin- 
ear elastic fields by incorporating short-scale physics 
through its additional dependence on the phase-field. 

We now consider a region SI that contains the process 
zone (crack tip) and write the integral of the divergence 
of the GEM tensor as a contour integral 
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We have decomposed the boundary of SI into: (i) a large 
loop {A B) around the tip in the unbroken material, 
where A {B) is at a height h below (above) the crack axis 
that is much larger than the process zone size but much 
smaller than the radius R of the contour, ^ <C ^ i?, 
and (ii) the segment {B A) that traverses the crack 
from B to A behind the tip, as illustrated in Fig. In 
both integrals, ds is the contour arclength clement and 
Uj the components of its outward normal. 

Eq. Q provides the basis to predict the crack speed 
and its path for quasi-static fracture. The Fi can be in- 
terpreted as the sum of all forces acting on the crack tip 
parallel {i = I) and perpendicular [i = 2) to the crack 
direction, including configurational, cohesive, and dissi- 
pative forces, represented by the three integrals terms 
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FIG. 1: Spatially diffuse crack tip region with (j> = 1/2 contour 
separating broken and unbroken material (thick solid line). 

from left to right, respectively. In the unbroken mate- 
rial where is constant, the tensor reduces identi- 
cally to the energy-momentum tensor introduced by Es- 
helby to compute the configurational force on the crack 
tip treated as a defect in a linear elastic field , and used 
in subsequent attempts to derive criteria for crack prop- 
agation and stability flHIlllll^ . Thus, the first integral 
in Eq. yields the configurational forces 

/ dsTijUj^G, / dsT2jnj ^Ge{0), (10) 

J A-*B Ja^B 

in the limit h/R 0. The first is the crack extension 
force and also Rice's J integral 13]. The second is the 
Eshelby torque Ge = dG{9)/de 9], where G{9) is the ex- 
tension force at the tip of a crack extended at a vanish- 
ingly small angle 9 from its local direction. This torque 
tends to turn the crack in a direction that maximizes G. 

An important new ingredient of the present derivation 
is the second portion of the line integral (/g_^^) of the 
GEM tensor that traverses the crack. This integral repre- 
sents physically the contributions of cohesive forces inside 
the process zone. To see this, we first note that the pro- 
files of the phase-field and the three components of the 
displacement can be made to depend only on X2 provided 
that the contour is chosen much larger than the process 
zone size and to traverse the crack perpendicularly from 
B to A. With this choice, we have that ni = —1, n2 — 0, 
along this contour and therefore that, for i = 1 

/ dsTijnj = - dx2Tii = -2j, (11) 

where the second equality follows from the fact that spa- 
tial gradients parallel to the crack direction (ditp^) give 
vanishingly small contributions in the limit /i/C — ^ +oo 
and i?/^ ^ -l-oo with h/R^ 0. This yields the expected 
result that cohesive forces exert a force opposite to the 
crack extension force with a magnitude equal to twice the 
surface energy. An analogous calculation for i — 2 yields, 
in the same limit ^ h ^ R, the other component of 
the force perpendicular to the crack direction 

/ dsT2jnj = -[ dx2T2i = ~2-ie{Q). (12) 

Jb-*A J-h 



This force is the direct analo g of the Herring torque = 
d'y/d9 on grain boundaries [lj|. This torque tends to 
turn the crack into a direction that minimizes the surface 
energy. Substituting the results of Eqs. (|l()|l to (|12|l into 
Eq. we obtain the two conditions 

Fi - G - - /i = 0, (13) 
F2 = Ge(0)-G,e(0)-/2 = 0, (14) 

where we have used the fact that Gee — '^le , and defined 
the dissipative forces 

/ + 00 P + OO 
j dxidx2 di(t)di(j), (15) 

by letting the area f2 tend to infinity since the integrand 
vanishes outside the process zone. Eq. ifT^ predicts the 
crack speed V ~ x(G — Gc)/ / dx{di4ioY foi' G close 
to Gc where (pQ is the phase-field profile for a station- 
ary crack 0. Eq. (O, in turn, predicts the crack path. 
G and Gg can be generally obtained from Eq. (|10|l . us- 
ing the known forms of the displacement fields near the 
tip. The J integral yields Eq. ((SJ for G. Ge can also be 
obtained directly from the expression for G{9). The lat- 
ter is instructive here to highlight important differences 
between plane strain and antiplane shear. Consider a 
straight crack parallel to the xi axis with stress intensity 
factors Ki and K2- Now extend this crack by a length L 
at a small angle 9 from this axis. The new stress in- 
tensity factors are given by Kl K-[ — "6X29/2 and 
K2 ~ K2 + Ki9/2 to linear order in 9 [l^] independent of 
L. Using Eq. Q with these new stress intensity factors 
to define G{9), we obtain at once Ge{0) = —2aKiK2- 
Substitution in Eq. H14|l . provides the condition 

K2 = -{G,em + .f2)/{2aKi), (16) 

which determines the crack path. In an isotropic mate- 
rial, this condition reduces to the principle of local sym- 
metry since Gee vanishes trivially, and f2 — 0- The 
latter follows from the symmetry of the inner phase- 
field solution for a propagating crack with K2 = 0, 
(f){xi,X2) ~ (j){xi, -'X2), which implies that the product 
di(j)d24> in Eq. (|15(l is anti-symmetric. In an anisotropic 
material, however, (jj is generally not symmetrical and /2 
only vanishes in the zero velocity limit where G ^ Gc- 

The same procedure can be repeated for pure an- 
tiplane shear where = if 3 — bnA2VL9 to linear or- 
der in 9 where 6 is a numerical constant, and hence 
Ge(0) '--^ 7^3^12 One important difference with plane 
loading is the divergence of Ge (0) with the crack exten- 
sion length L. This divergence is also reflected in a 
dependence of the integral in Eq. H10|l on the radius R of 
the contour enclosing the tip. Since the only natural cut 
off for this divergence is the system size, this result seems 
to imply that the crack path cannot be predicted solely 
in terms of local conditions at the tip for mode III . We 
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expect, however, this divergence to be cut off in a real ex- 
periment by the process zone size ^ due to the irreversible 
nature of the fracture process. Namely, fracture surfaces 
at a distance behind the tip larger than ^ should be essen- 
tially immobile, which implies that Ge{0) ~ -f^3^2-\/C up 
to a numerical prefactor. Eq. H14I) then yields the local 
symmetry condition A2 — in the isotropic limit, where 
the symmetry of the phase-field profile makes /2 vanish, 
as explained above. We will present elsewhere numerical 
results that validate this condition for mode III. We focus 
in the remainder of this letter on plane strain. 

We use a simple anisotropic extension of the phase-field 
model of Ref. '61 with an energy density 

(|V0|2 + edi(j)d2^) /2 + g{(j)) {£ 

strain £c), (17) 

where u.^ = (diUj + djUi)/2 is the strain tensor and 
^^stiain = /iufj is the strain energy. No asymmetry 

between dilation and compression is included since this 
is not necessary here to test our predictions. The broken 
state of the material becomes energetically favored when 



^strain cxcceds a threshold 8c and g{(t)) — 
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a monotonously increasing function of that describes 
the softening of the elastic energy at large strain. By 
repeating the analysis of Ref. 0], we obtain that 7(0) = 
— (e/2) sin 20 where 7 reduces to the isotropic sur- 
face energy of Ref. in the e ^ limit. 

We test our prediction for the initial angle of a kink 
crack. Eq. is solved numerically using an Euler ex- 
plicit scheme to integrate the phase-field evolution and 
a successive over relaxation (SOR) method to calculate 
the quasi-static displacement fields Ui and U2 at each 
time step. We used as initial condition a straight hori- 
zontal crack of length 2W centered in a strip of length 
AW horizontally and 2W vertically, with fixed values of 
Ml and U2 on the strip boundaries that correspond to 
the singular stress fields defined by Eq. for prescribed 
values of Ki and K2- We used A/^ ^ \ [a — 3/(8^)], a 
grid spacing Axi = Aa;2 = 0.1^, and W — 50^, where 
the process zone size ^ = ^fnjJJiE^. We checked that 
the results are independent of width and grid spacing. 

We have verified that the kink angle is well pre- 
dicted by the local symmetry condition K2 = in the 
isotropic limit, which implies that d ~ —2K2/K1. In the 
anisotropic case, we choose K2 = G and G just slightly 
above Gc such that /2 can be neglected in Eq. . Sub- 
stituting K2 ~ K\d 12 in Eq. (|16|l and using the fact that 
(1 - v'^)Kl/E « 27(0) for G close to G^ we obtain the 
prediction for the kink angle 9 ~ —70/7 ~ e/2 which 
is strictly valid for <C 1 and G ^ Gc- This predic- 
tion is in good quantitative agreement with the results of 
phase-field simulations as shown in Fig. |2 

An interesting implication of our results for crystalline 
materials is that \K2\ should exceed some threshold for 
a cleavage crack to change direction. Using the expected 
cusp behavior of the surface energy for small angle near a 
cleavage plane, 7(6') = 7o(l-l-(5|6'| -I- . . .), Eq. predicts 




FIG. 2: Kink angle 6 versus surface energy anisotropy e pre- 
dicted as 6 = e/2 and extracted from phase-field simulations 
(filled circles) for G/Gc ~ 1.1. Inset: phase-field simulation 
for e — 1.2 (0 = 1/2 contours are equally spaced in time). 



that this threshold is E^qS/[{1 — l''^)Ki] for G ~ Gc since 
this equation cannot be satisfied for any smaller value of 
K2 for small. It should be hopefully possible to test 
this prediction experimentally as well as to explore the 
validity of this new condition on K2 for curvilinear paths. 
The extension of the present analysis to include inertial 
effects and to three dimensions where fracture paths are 
geometrically more complex is an important future direc- 
tion. Work along this line is presently in progress. 
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